ACP focus stress biotic
Solène Pety
25/03/2021
library(FactoMineR)
library(factoextra)
library("corrplot")
library(ggsci)
library(plotly)
library(ggplot2)
library(dplyr)
library(tidyverse)
library(ComplexUpset)
library(pander)
library(knitr)
library(discoveR)
library(ade4)
library(adegraphics)
library(vegan)
library(MASS)
library(ComplexUpset)
library(RColorBrewer)Dans tout le jeu de données, on retrouve :
- 9 stress biotiques (Biotrophic bacteria, Fungi, Nématodes, Oomycète, Rhodococcus, Stifenia, Necrotrophic bacteria, Virus et Other biotic).
- 9 stress abiotiques (Heavy metal, UV, Drought, Gamma, Nitrogen, Oxydative stress, Salt, Temperature et Other abiotic).
Seuls les stress biotiques sont étudiés ici.
Les fichiers utilisés ont deux colonnes informatives avec le stress appliqué sur l’échantillon en première colonne et le SWAP_ID pour pouvoir remonter aux informations de l’expérience précise.
Le fichier de chaque Gene Set est chargé à partir de la liste fournie dans le script. Les échantillons pour ce set sont séparés enuite pour les stress biotiques et abiotiques.
# Chemin vers dossier fichiers .txt et chargement
setwd("../docs_txt/")
data=read.table('Gene_Swap_NO_NA.dat',header=TRUE,sep='\t')
# Noms des colonnes sans les variables qualitatives
Gene=noquote(colnames(data[,c(-1,-2,-3)]))
#Noms des fichiers pour lecture dans boucle
GOSLIM=c("circadian_rythm.txt","abiotic.txt","biotic.txt","endogenous.txt",
"external.txt","flower.txt","growth.txt","light.txt","photo.txt",
"stress.txt")
matrix=cbind(rep(0,17341),Gene) #initialisation matrice
#Remplacer 0 par TRUE ou FAlSE pour chaque gene set
for(i in 1:10){
datago=read.table(GOSLIM[i],header=TRUE,sep='\t')
colgo=noquote(colnames(datago[,c(-1,-2,-3)]))
test=match(Gene,colgo,nomatch =0)
test[test!=0]<-"TRUE"
test[test==0]<-"FALSE"
matrix=cbind(matrix,test)
}
matrix=as.data.frame(matrix)
matrix=matrix[,-1]#enlever la première colonne de zéro d'initialisation
names(matrix)=c("Gene","Circadian","Abiotic","Biotic","Endogenous_stimulus",
"External_Stimulus","Flower","Growth","Light","Photosynthesis",
"Stress")
goslim=colnames(matrix)[2:11]
#Construction upset
(upset(matrix,goslim, name="GO SLIM", min_size=25,min_degree=1)
+ ggtitle("Tous les gènes (173441)"))
setwd("../docs_txt/")
#Variables importantes pour code, noms des fichiers pour lecture, nom des gene set pour entêtes personnalisées, différents stress..
files=c("circadian_rythm.txt","abiotic.txt","biotic.txt","endogenous.txt",
"external.txt","flower.txt","growth.txt","light.txt","photo.txt",
"stress.txt","random.txt")
name=c("Circadian rythm","Abiotic","Biotic","Endogenous stimulus",
"External stimulus","Flower","Growth","Light","Photosynthesis",
"Stress","Random")
biotic=c("BIOTROPHIC.BACTERIA","FUNGI","NECROTROPHIC.BACTERIA","NEMATODES","OOMYCETE","OTHER-BIOTIC","RHODOCOCCUS","STIFENIA","VIRUS")
#Boucle principale lecture Gene set
for(i in 1:11){
#Chargement Gene set
data=read.table(files[i],header=TRUE,sep='\t')
#Conservation des échantillons seulement abiotiques
data_biotic=data[data[,2] %in% biotic, ]
row.names(data_biotic) <- data_biotic$Swap_ID
#ACP avec données abiotiques
res1=PCA(data_biotic[,c(-1,-3)],
scale.unit=FALSE,
graph=FALSE,
quali.sup=1)
fic.pca1<-dudi.pca(data_biotic[,c(-1,-2,-3)],
center=FALSE,
scale=FALSE,
scannf=FALSE)
#BCA données abiotiques
fic.bca1<-bca(fic.pca1,fac=as.factor(data_biotic$vec),scannf=FALSE)
fic.tst<-randtest(fic.bca1)
#--------------------------------------------------------------------------------------------------#
#Ecriture dans ouput HTML
#--------------------------------------------------------------------------------------------------#
cat("\n")
pander::pandoc.header(paste0("**Gene set : ",name[i]," (",dim(data)[2]-2," gènes)**\n"), level = 1)
#--------------------------------------------------------------------------------------------------#
#ACP
#--------------------------------------------------------------------------------------------------#
cat("\n")
pander::pandoc.header(paste0("**ACP (SWAP_ID)**\n"), level = 2)
g1<-fviz_pca_ind(res1,
repel=TRUE,
habillage = 1, # color by groups
addEllipses = TRUE, # Concentration ellipses
palette="Set1")
print(g1)
cat("\n")#Screeplot intertie
pander::pandoc.header(paste0("*Inertie*\n"), level = 3)
print(fviz_screeplot(res1, addlabels = TRUE, ylim = c(0, 50)))
cat("\n")#Screeplot contribution individus et variables
pander::pandoc.header(paste0("*Contribution*\n"), level =3)
cat("\n")
pander::pandoc.header(paste0("**Individus**\n"), level =4)
print((fviz_contrib(res1, choice = "ind", axes = 1, top = 10)
+
fviz_contrib(res1, choice = "ind", axes = 2, top = 10)))
cat("\n")
pander::pandoc.header(paste0("**Variables**\n"), level =4)
print((fviz_contrib(res1, choice = "var", axes = 1, top = 10)
+
fviz_contrib(res1, choice = "var", axes = 2, top = 10)))
cat("\n")#Boxplot des meilleurs contributeurs
pander::pandoc.header(paste0("**Distribution meilleurs contributeurs**\n"), level =4)
#Récupération des labels des contributeurs via les plots factoextra
#Découpage en 2 fois pour les deux dimensions
p=fviz_contrib(res1, choice = "var", axes = 1, top = 10)
p2=fviz_contrib(res1, choice = "var", axes = 2, top = 10)
contrib=p$data[order(p$data$contrib,decreasing=TRUE), ]
contrib2=p2$data[order(p2$data$contrib,decreasing=TRUE), ]
contrib=contrib[1:5,]
contrib2=contrib2[1:5,]
#Match entre noms des colonnes data_abiotic et les noms des meilleurs contributeurs
col_names=colnames(data_biotic)
index=match(contrib$name,col_names)
index=index-2
index2=match(contrib2$name,col_names)
index2=index2-2
#Récupération des données de puces (sans colonnes qualitatives) pour les colonnes d'intérêt
contrib_data=data_biotic[,c(-1,-2)]
contrib_data=contrib_data[,index]
contrib_data2=data_biotic[,c(-1,-2)]
contrib_data2=contrib_data2[,index2]
data_long=contrib_data%>%pivot_longer(
cols = c(1:5),
names_to = "Gene",
values_to = "Expression")
g1<-ggplot(data_long)+
geom_boxplot(aes(x=Gene, y = Expression))+
theme(axis.text.x = element_text(angle = 60, hjust = 1))+
ggtitle("Contributions Dim 1")
data_long2=contrib_data2%>%pivot_longer(
cols = c(1:5),
names_to = "Gene",
values_to = "Expression")
g2<-ggplot(data_long2)+
geom_boxplot(aes(x=Gene, y = Expression))+
theme(axis.text.x = element_text(angle = 60, hjust = 1))+
ggtitle("Contributions Dim 2")
#Détermination des limites en fonction des valeurs des deux dimensions pour rendre comparaison plus facile
g1<-g1+ylim(min(data_long$Expression,data_long2$Expression),
max(data_long$Expression,data_long2$Expression))
g2<-g2+ylim(min(data_long$Expression,data_long2$Expression),
max(data_long$Expression,data_long2$Expression))
print(g1+g2)
#--------------------------------------------------------------------------------------------------#
#BCA
#--------------------------------------------------------------------------------------------------#
cat("\n")
pander::pandoc.header(paste0("**BCA**\n"), level = 2)
cat("\n")
pander::pandoc.header(paste0("*Toutes positions*\n"), level = 3)
cat("\n")
.z<-paste("BCA sur position \n ratio= ",round(fic.bca1$ratio,2)," pval= ",signif(fic.tst$pvalue,digits=2))
s.class(fic.bca1$ls,
fac=as.factor(data_biotic$vec),
col=rainbow(9),
psub=list(text=.z,cex=1,position="topleft"),
plabels.cex=0,
key.cex=0.8,
key.text.cex=0.8)
cat("\n")
pander::pandoc.header(paste0("*Par position*\n"), level = 3)
#Plot pour chaque stress
plotlist=list()
for(j in 1:9){
s.class(fic.bca1$ls,
fac=as.factor(data_biotic$vec),
psub=list(text=paste("BCA \n", biotic[j]),cex=1,position="topleft"),
plabels.cex=0,
key.cex=0.4,
key.text.cex=0.4,
col=c(rep(grey(0.95),j-1),brewer.pal(n=9,name="Set1")[j],rep(grey(0.95),9-j)),
plot=FALSE,
chullSize=0,
ellipseSize=0)->plotlist[j]
}
ADEgS(plotlist,layout=c(3,3))
}Gene set : Circadian rythm (58 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : Abiotic (662 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : Biotic (421 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : Endogenous stimulus (524 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : External stimulus (553 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : Flower (187 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : Growth (223 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : Light (293 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : Photosynthesis (76 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : Stress (1178 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
Gene set : Random (51 gènes)
ACP (SWAP_ID)
Inertie
Contribution
Individus
Variables
Distribution meilleurs contributeurs
BCA
Toutes positions
Par position
